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1. Introduction 


Modeling the inviscid air flow and its constituents over a hypersonically flying body 
requires a large system of Euler and chemical rate equations in three spatial coordinates. In 
most cases, the simplest approach to solve for the variables would be based on explicit inte- 
gration of the governing equations. But the standard techniques are not suitable for this 
purpose because the integration step size must be inordinately small in order to maintain 
numerical stability. The difficulty is due to the stiff character of the difference equations, as 
there exists a large spectrum of spatial and temporal scales in the approximation of physical 
phenomena by numerical methods. For instance, in the calculation of gradients caused by 
shock and by cooled wall on a coarse grid, unchecked numerical errors eventually will lead to 
violent instability, and in calculations of species near chemical equilibrium, a small error in 
one species will give rise to a large error in the source term for other species. Despite the 
different nature of the stiffness in a complex system of equations, the most effective approach 
is believed to be implicit integration. The step increment is no longer dictated by the stabil- 
ity criteria for explicit methods, but instead is dictated by the degree of linearization 
introduced to the governing equations and by the order of desired accuracy. The lineariza- 
tion is enacted by means of Jacobian matrices, resulting from the differentiation of the flux 
as well as the rate production terms with respect to dependent variables. The backward 
Euler scheme is then applied to discretize the partial differential equations and to convert 
them into a system of linear difference equations in vector form. As this particular approach 
has the A-stable property, it is the one recommended by Lomax and Baileyd) for one- 
dimensional nonequilibrium flow studies. However, in the practice of solving flow problems 
in multidimensions, it was not clear then how to deal with the mammoth size of the sparse 
block matrix equations. The implementation of an implicit method in the solution procedure 
could be as prohibitively expensive as a modified Runge-Kutta method.(2) 


2. Formulation 


In view of the drawbacks associated with the implicit methods, other concepts have 
been evolved to avoid using the fully coupled approach. The most notable concepts probably 
are the hybrid explicit-implicit techniques that focus on the minimization of the stiffness due 
to the chemical production(3-4) and the time-split explicit method devised in such a manner 
that the flow and species equations are integrated separately according to the stability 
criteria of each.(5) These plausible ideas have provided limited success in two-dimensional 
problems, for which the flow equations are not stiff. The dilemma in regard to the unre- 
stricted stability conditions and the excessive manipulation of a matrix- vector equation may 
be resolved now by taking a different path in seeking an efficient method. 

The objective of the research reported herein is to evaluate the efficiency and robust- 
ness of two variants of an implicit method that has recently been developed to investigate 
the complete flowfield around an entry vehicle. The most challenging aspect of nonreacting 
flow computation was the generation of a computational grid so as to allow a single data set 
structure for both shock layer and trailing wake in one zone. The Euler equations were cast 
in the domain between the wall and the outer surface consisting of the bow shock and were 
solved by a modified version of the ADI factorization technique.(6) The methodology can be 
extended to consider reacting species in several ways; namely, coupled flow and species, 
decoupled flow from species and simultaneous species solution, and decoupled successive 
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species solution. The differences among the three approaches are mainly in the degree of 
linearization and the ensuing amount of computation. The coupled approach would be the 
favorable one to use if the initial conditions are close to the final solution and the temporal 
accuracy is critical. The last two approaches are more appealing whenever the final solution 
is obtained after a reasonable, economical number of iterations. In the hypersonic flowfield 
with large subsonic regions surrounding the body, the initialization is at best very crude. 
Hence, the issue of efficiency is not as important as the robustness issue. Besides, the third 
approach has the potential of being the least computationally intensive method as the 
ionization is considered (11 vs 7 species). 


3. Outline of the Decoupled Methods 

The ideas can be elucidated by the one-dimensional problem having two species. 
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and A is the duct area; the subscripts M and A refer to molecular and atomic species respec- 
tively; and co is the production term. Standard notation is used otherwise. 

An alternate equation equivalent to Eq. (1) is 
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It is readily seen from Eq. (2) that a weak relationship exists between (p,u, e) and (C M ,C A ). 
The loosely coupled property in P and A will be exploited further in the following section. 
Since neither Eq. (1) nor Eq. (2) is sufficient to describe the flowfield, two equations of T and 
P 





(3) 
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also are needed. Relations of the conservation of mass concentration, Cf, and of charges are 
used for more complicated chemistry. 

The production term, to^ = oip(g,T,C£), is proportional to (C f - C^ e )/x, where Co e is the 
equilibrium value and x is the reaction rate. It can have astronomical value even when nor- 
malized by the flow resident time in some portion of the flow. To cope with such stiffness, its 
relationship with the dependent variable should be analyzed and included in the algorithm 
by means oi a Taylor series expansion as follows: 
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The superscript ° denotes the known values, and A represents the difference between the 
present and the known values. By incorporating Eq. (4) to Euler’s backward scheme, Eq. (2) 
is converted to difference form in shorthand notation. 


M k AV* +1 + A<8 (A* AV* +1 ) = ffffS 

l l X \ l l ] 

where 8 X refers to the centered-difference operator, i is the index of the grid and 


RHS = -A tP~ x 



(5) 
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The appearance of M in Eq. (5) and of the partial derivatives in the elements of M 
serves the role of tempering the resultant A V t . Simplifications can be made to M by neglect- 
ing off-diagonal elements or only those elements associated with ( )p and ( )e. Then, Eq. (5) is 
effectively decoupled into two groups for flow and species variables. 

M has two splitup versions as follows: 



It is interesting to note that, in the simplest version, Me' 1 multiplies the RHS and the 
operator 8^ and results in the reduction of the time increment Af by a factor of 1 + l(co • ) Ci AH. 
Since each species equation is then integrated by a different At, this approach uses the same 
idea advocated in Li.C5) 

Equation (5) is a tridiagonal block or scalar system of equations depending on the form 
of M for which the solution procedure is selected. It represents one of the three steps in 
solving three-dimensional problems. The truncation error of the method is (Af , Ax 4 * 6 ). By 
contrast, higher order explicit methods can be constructed by using 

A V k+1 = M~ l (RHS) (6) 

I l 

although At may be restricted to a narrow corridor of AA t, where A is the largest eigenvalue 
of A. 


4. Discussion of Results 

The development of decoupled implicit methods was based on two body configurations 
of practical interest and under flow conditions that give rise to extreme levels of dissociation 

and ionization. The fully implicit method has not been implemented in a computer code 
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because it appears to be prohibitively costly for an 11-species model on today’s scalar 
computer. The simultaneous species solver was used primarily to compare its performance 
with the less expansive successive species solver. 

The first case was a sphere of R N = 0.328 ft, at M = 10 and h = 100 kft. Some results 
of M = 10 reacting flow over the frontal portion of the sphere (table 3.16 in BelotserkoskiyU)) 
were used to verify the accuracy of present results. In as much as the body is simple and the 
dissociation is weak, a coarse grid of 6 X 20 was found to be adequate. 

Figure 1 shows the convergence history of the mass fraction of the oxygen molecule in 
terms of the maximum incremental and the stagnation values. As shown, both figures are 
needed to exhibit the realistic rate of convergence. The noticeably slow relaxation process in 
chemical nonequilibrium is indicative of the mutual interactions between the convection and 
the production of species allowed by the step size, CN = 2. Attempts to increase the Courant 
number to CN = 4 failed to control the growth of error; indeed, a wide band of fluctuation 
was observed in AC 02 - After ACO 2 reac hed and stayed at the level of truncation error, 

CO 2 began to converge at the forward stagnation point. Using the convergence history 
obtained from both methods, it was estimated that the successive method is a factor of two 
slower. 

Figure 2 presents the temperature and species distributions across the shock layer in 
three angular orientations: 6 = 0, n/2, and n. An excellent agreement is found in tempera- 
ture, but discrepancies are seen in the species distributions. As pointed out in Li, (5) the rate 
constants usually exert stronger influences on the species than on the temperature. 

The second case investigated was an aerobrake of R^ = 20 ft at M = 34.8 and h = 250 
kft. The computational grid, 13 X 36, was similar to that developed in Li(6) for nonreacting 
flow computations. A comparison of the convergence history between the simultaneous and 
the successive methods is shown in Fig. 3 for the nitrogen atom. The formal method has a 
faster rate in the first 100 iterations, but commands as many iterations as the latter method 
to reach the asymptotic value. The simultaneous method has produced a great deal of 
fluctuation that is atypical for flow with high levels of dissociation. 

The successive method was slightly more efficient (by 15%) than the simultaneous 
method. Out of the total execution time, about 50% was spent to solve for the flow variables, 
43% for the 5X5 (Hq matrix, and the rest of the execution time for the species. During each 
iteration, the analytical procedure determining co and a c was called twice for each grid point 
such that no extra core was assigned to store them. 

The temperature contours obtained for the two cases are exhibited in Fig. 4. A com- 
plete flowfield can be obtained with 200 to 400 iterations on the coarse grid laying between 
the wall and the outer surface. 

Figure 5 illustrates the convergence history of CV at the front and the rear stagnation 
points. This is the same aerobrake considered in case 2, but at a 20° angle of attack. The 
nonequilibrium flow in the shock layer changed little after 200 iterations, yet, in the near 
wake, probably more than 400 iterations were required to ensure convergence. Note that the 
stagnation value of Cy was lower than that predicted for the zero angle of attack case. 

Figure 6 gives the contour plot of Mach and e~ (no./cm 3 ) made on the pitch plane. The shear 
layer was not clearly visible because the grid of 13 X 36 X 6 was quite coarse. The distribu- 
tion of electrons displays three orders of magnitude difference between the shock layer and 
the wake. 


5. Conclusion 

It was found that the conventional implicit techniques can effectively reduce the stiff- 
ness associated with the chemical production term in the rate equations. The successive 
solution for the species was as stable as the simultaneous solution. The fact that reactive 
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Jacobians were used to scale down the time increment was essentially regarded as a means 
to separate the rate equations and to perform integrations according to individual reaction 
time. In summary, the numerical procedure consists of decoupling the flow from species 
equations, solving them by block and scalar tridiagonal procedures, respectively, and 
utilizing the factorization ADI technique to tackle two- or three-dimensional problems. The 
procedure is stable and the computation time of species varies linearly with the number of 
species. Significant reduction of computation time can be achieved by storing the reactive 
Jacobians and by updating them periodically in the coarse of iterations. 
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Fig. 1 Comparison of the convergence history between the simultaneous and successive 
solutions for case 1: (a) maximum value of the incremental Cc> 2 > (b) forward stagnation 
value of Co 2 > normalized by its maximum value. 



(a) v (b) 


Fig. 2 Comparison of temperature and species mass fractions: (a) temperature 
distributions vs normalized distance from outer surface to wall; (b) species distributions 
along the forward stagnation line between shock and wall. 
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Fig. 3 Comparison of the convergence history between the simultaneous and successive 
solutions for case 2: (a) maximum value of the incremental C N ; (b) forward stagnation 
value of C N , normalized by its normalized value. 




Fig. 4 Temperature contours for cases 1 and 2: (a) casel;(b) case 2. 



Fig. 5 Convergence history of case 2 at a = 20°: (a) front stagnation value of C^; (b) rear 
stagnation value of C^. Both are normalized quantities. 


8 





Fig. 6 Contours of Mach number and electron number density (dash line = sonic line): (a) 
Mach number; (b) e- (no./cm3). 
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